

## #######################################################################################################
## 
## Replication R command file for:
## Graeme Blair and Kosuke Imai, "Statistical Analysis of List Experiments," Political Analysis, In Press.
##
## Created 27 October 2011
## 
## Contact Graeme Blair <gblair@princeton.edu> with questions
## 
## #######################################################################################################

## ###
## NB: the code here uses functions of the R list package available through CRAN at
## http://cran.r-project.org/web/packages/list/index.html
## Version 4.1 of the software is included in this replication archive,
## which was used to confirm that this code replicates the tables and figures in the paper.
## ###

## Figure 1 - "Differences between the Binomial and Poisson-Binomial Distributions"
## pg 15

rm(list=ls())

require(list)

dpoisbinom <- function(y, p) {
				
  k <- length(p)
  return(.C("dpoisbinom", as.integer(y), as.integer(k), as.double(p), res = double(1), PACKAGE = "list")$res)
  
}

dpoisbinomE <- function(y, p) {

  return(R.poisbinomE(y, p)*prod(1-p))
  
}

R.poisbinomE <- function(k, p) {

  maxk <- max(k)
  return(.C("RpoisbinomEff", as.integer(maxk), 
            as.double(p), as.integer(length(p)), Rs = double(maxk+1),
            PACKAGE = "list")$Rs[k+1])

}

p1 <- c(.9, .7, .5, .3, .1)
p2 <- c(.7, .6, .5, .4, .3)
p3 <- c(.8, .5, .45, .4, .35)
p4 <- c(0.5, 0.4, 0.3, 0.2, 0.1)
p5 <- c(0.75, 0.3, 0.2, 0.15, .1)

x <- 0:5

binomdens5 <- dbinom(x, size = 5, p = .5)
binomdens3 <- dbinom(x, size = 5, p = .3)
dens1 <- dpoisbinomE(x, p1)
dens2 <- dpoisbinomE(x, p2)
dens3 <- dpoisbinomE(x, p3)
dens4 <- dpoisbinomE(x, p4)
dens5 <- dpoisbinomE(x, p5)

binom.cum5 <- rep(0, length(binomdens5))
binom.cum3 <- rep(0, length(binomdens3))
d1.cum <- rep(0, length(p1))
d2.cum <- rep(0, length(p2))
d3.cum <- rep(0, length(p3))
d4.cum <- rep(0, length(p4))
d5.cum <- rep(0, length(p5))
for(i in 1:length(dens1)){
	if(i>1) last <- i-1
	if(i==1) last <- i
	binom.cum5[i] <- binom.cum5[last] + binomdens5[i]
	binom.cum3[i] <- binom.cum3[last] + binomdens3[i]

	d1.cum[i] <- d1.cum[last] + dens1[i]
	d2.cum[i] <- d2.cum[last] + dens2[i]
	d3.cum[i] <- d3.cum[last] + dens3[i]
	
	d4.cum[i] <- d4.cum[last] + dens4[i]
	d5.cum[i] <- d5.cum[last] + dens5[i]

}

## Generate the plots

pdf("figure1.pdf", width = 8, height = 4.5, paper = "special")
par(mfcol = c(1, 2), cex = .8) #, mar = c(2, 3, 2, 1), oma = c(2.5, 2.5, 2.5, 0))

bp <- plot(c(0:5), dens4, type = "n", xlim=c(-.25, 5.25), ylim = c(0,.5), axes = F, xlab = "Count", ylab = "Density", main = "mean probability = 0.5", omi = 1)

for(i in 0:(length(x)-1)) rect((i-0)-.42,0, (i-0)+.42, binomdens5[i+1], dens = 0, border = "gray35", lwd = 1.5)
for(i in 0:(length(x)-1)) rect((i-.2)-.13,0, (i-.2)+.13, dens1[i+1], col = "gray30", border = NA)
for(i in 0:(length(x)-1)) rect((i+.2)-.13,0, (i+.2)+.13, dens2[i+1], col = "gray60", border = NA)

axis(1)
axis(2)

legend(-.1, .53, c("Binomial with p = 0.5","Poisson-Binomial with p = (0.9, 0.7, 0.5, 0.3, 0.1)","Poisson-Binomial with p = (0.7, 0.6, 0.5, 0.4, 0.3)"), fill = c(NA, "gray30", "gray60"), cex = .8, border = c("gray35", "gray30", "gray60"), bty = "n")

bp <- plot(c(0:5), dens4, type = "n", xlim=c(-.25, 5.25), ylim = c(0,.5),  axes = F, xlab = "Count", ylab = "Density", main = "mean probability = 0.3")

for(i in 0:(length(x)-1)) rect((i-0)-.42,0, (i-0)+.42, binomdens3[i+1], dens = 0, border = "gray35", lwd = 1.5)
for(i in 0:(length(x)-1)) rect((i+.2)-.13,0, (i+.2)+.13, dens4[i+1], col = "gray60", border = NA)
for(i in 0:(length(x)-1)) rect((i-.2)-.13,0, (i-.2)+.13, dens5[i+1], col = "gray30", border = NA)

axis(1)
axis(2)

legend(-.1, .53, c("Binomial with p = 0.3","Poisson-Binomial with p = (0.75, 0.3, 0.2, 0.15, 0.1)","Poisson-Binomial with p = (0.5, 0.4, 0.3, 0.2, 0.1)"), fill = c(NA, "gray30", "gray60"), cex = .8, border = c("gray35", "gray30", "gray60"), bty = "n")

dev.off()


## ####
## Table 2 - "Observed Data from the List Experiment in the 1991 National Race and Politics Survey"
## pg 16

rm(list=ls())
require(list)

## data is included in the R list package loaded above
data(multi)

table(multi$y, multi$treat)


## ####
## Table 3 - "Estimated Coefficients from the Combined Logistic Regression Model ..."
## pg 17

rm(list=ls())
require(list)

# Fit EM algorithm ML model for multiple sensitive items

multi.results <- ictreg(y ~ male + college + age + south + south:age, treat = "treat", 
	      	 	J = 3, data = multi, method = "ml", 
			multi.condition = "level")

summary(multi.results)


## ####
## Figure 2 - "Estimated Proportions of Respondents Answering Each Sensitive Item in the Affirmative by Respondent Age"
## pg 18

rm(list=ls())
require(MASS)
require(list)

J <- 3

results <- ictreg(y ~ male + college + age + south + age:south, treat = "treat", 
	      	 	J = 3, data = multi, method = "ml", 
			multi.condition = "level")

par.both <- coef(results)

n.draws <- 10000

logistic <- function(x) exp(x)/(1+exp(x))

## set up data

X.fam <- with(race, cbind(1, male, college, age, south, age*south))
X.affirm <- with(affirm, cbind(1, male, college, age, south, age*south))

X.fam.south <- with(subset(race, south == 1), cbind(1, male, college, age, south, age*south))
X.affirm.south <- with(subset(affirm, south == 1), cbind(1, male, college, age, south, age*south))

X.fam.north <- with(subset(race, south == 0), cbind(1, male, college, age, south, age*south))
X.affirm.north <- with(subset(affirm, south == 0), cbind(1, male, college, age, south, age*south))

k <- 6

south.fam.multi <- north.fam.multi <-
  south.affirm.multi <- north.affirm.multi <- rep(NA, 88-18)
 
for (age in 18:88) {
  
  X.south.fam <- X.fam.south

  X.south.fam[,4] <- age/10
  X.south.fam[,6] <- X.south.fam[,5] * X.south.fam[,4]
  
  X.north.fam <- X.fam.north

  X.north.fam[,4] <- age/10
  X.north.fam[,6] <- X.north.fam[,5] * X.north.fam[,4]
  
  X.south.affirm <- X.affirm.south

  X.south.affirm[,4] <- age/10
  X.south.affirm[,6] <- X.south.affirm[,5] * X.south.affirm[,4]
  
  X.north.affirm <- X.affirm.north

  X.north.affirm[,4] <- age/10
  X.north.affirm[,6] <- X.north.affirm[,5] * X.north.affirm[,4]
  
  par.h <- par.both[(2*k+3):length(par.both)]
  par.g.fam <- par.both[(k+2):((k)*2+2)]
  par.g.affirm <- par.both[1:(k+1)]
  
  pred.north.fam <- matrix(NA, nrow = nrow(X.north.fam), ncol = J+1)
  pred.south.fam <- matrix(NA, nrow = nrow(X.south.fam), ncol = J+1)
  pred.north.affirm <- matrix(NA, nrow = nrow(X.north.affirm), ncol = J+1)
  pred.south.affirm <- matrix(NA, nrow = nrow(X.south.affirm), ncol = J+1)
  
  for (y in 0:J) {
    
    pred.north.fam[,y+1] <- logistic(cbind(X.north.fam,y) %*% par.g.fam) *
      dbinom(x = y, size = J, prob = logistic(X.north.fam %*% par.h), log = FALSE)
    pred.south.fam[,y+1] <- logistic(cbind(X.south.fam,y) %*% par.g.fam) *
      dbinom(x = y, size = J, prob = logistic(X.south.fam %*% par.h), log = FALSE)
    
    pred.north.affirm[,y+1] <- logistic(cbind(X.north.affirm,y) %*% par.g.affirm) *
      dbinom(x = y, size = J, prob = logistic(X.north.affirm %*% par.h), log = FALSE)
    pred.south.affirm[,y+1] <- logistic(cbind(X.south.affirm,y) %*% par.g.affirm) *
      dbinom(x = y, size = J, prob = logistic(X.south.affirm %*% par.h), log = FALSE)
    
  }
  
  south.fam.multi[age-17] <- mean(apply(pred.south.fam, 1, sum))
  north.fam.multi[age-17] <- mean(apply(pred.north.fam, 1, sum))
  
  south.affirm.multi[age-17] <- mean(apply(pred.south.affirm, 1, sum))
  north.affirm.multi[age-17] <- mean(apply(pred.north.affirm, 1, sum))
  
}

pdf("figure2.pdf", width = 6.5, height = 3)
par(mfrow = c(1,2), mar = c(4, 4, 1, 1), las = 1, oma = c(1.5, 0, 0, 0), cex = 0.6)
plot(18:88, south.fam.multi, type = "l", lty = 2, lwd = 1, ylim = c(0,1),
     ylab = "Estimated Proportion", xlab = "Age", axes = F, main = "Black Family")
lines(18:88, north.fam.multi, lty = 1, lwd = 1)
axis(side = 2, at = seq(from = 0, to = 1, by = 0.2))
axis(side = 1, at = seq(from = 20, to = 90, by = 10))
text( 22, 0.25, "South", srt = 10)
text( 26, 0.07, "Non-South", srt = 2)

plot(18:88, south.affirm.multi, type = "l", lty = 2, lwd = 1, ylim = c(0,1), axes = F,
     ylab = "Estimated Proportion", xlab = "Age", main = "Affirmative Action")
lines(18:88, north.affirm.multi, lty = 1, lwd = 1 )
axis(side = 2, at = seq(from = 0, to = 1, by = 0.2))
axis(side = 1, at = seq(from = 20, to = 90, by = 10))
text( 22, 0.85, "South", srt = -9)
text( 25, 0.23, "Non-South", srt = 25)

dev.off()


## ####
## Figure 3 - "Estimated Proportions of Respondents Answering the Sensitive Item in the Affirmative by Partisanship"
## pg 18

rm(list=ls())
require(MASS)
require(list)

mis.list <- subset(mis, list.data == 1 &
                     (democrat==1 | republican==1 | independent==1))

mis.sens <- subset(mis, sens.data == 1 &
                     (democrat==1 | republican==1 | independent==1))

fit.list <- ictreg(y ~ age + college + male + south + democrat + republican,
                   J = 4, data = mis.list, method = "ml", fit.start = "lm")

fit.sens <- glm(sensitive ~ age + college + male + south + democrat +
                     republican, data = mis.sens, family = binomial("logit"))

vcov.list <- vcov(fit.list)
par.list <- coef(fit.list)

vcov.sens <- vcov(fit.sens)
par.sens <- coef(fit.sens)
              
n.draws <- 10000

logistic <- function(x) exp(x)/(1+exp(x))

X.list <- as.data.frame(with(mis.list,
                             cbind(1, age, college, male, south, democrat, republican)))

X.sens <- as.data.frame(with(mis.sens,
                             cbind(1, age, college, male, south, democrat, republican)))

k <- ncol(X.list)

X.democrat.list <- as.matrix(subset(X.list, democrat == 1))
X.republican.list <- as.matrix(subset(X.list, republican == 1))
X.independent.list <- as.matrix(subset(X.list, democrat == 0 & republican == 0))

X.democrat.sens <- as.matrix(subset(X.sens, democrat == 1))
X.republican.sens <- as.matrix(subset(X.sens, republican == 1))
X.independent.sens <- as.matrix(subset(X.sens, democrat == 0 & republican == 0))

democrat.list <- republican.list <- independent.list <-
  democrat.sens <- republican.sens <- independent.sens <-
  democrat.diff <- republican.diff <- independent.diff <- rep(NA, n.draws)

draws.list <- mvrnorm(n = n.draws, mu = par.list[1:k], Sigma = vcov.list[1:k,1:k])
draws.sens <- mvrnorm(n = n.draws, mu = par.sens, Sigma = vcov.sens)

for (d in 1:n.draws) {
  
  par.g <- draws.list[d, ]
    
  pred.democrat.list <- logistic(X.democrat.list %*% par.g)

  pred.democrat.sens <- logistic(X.democrat.sens %*% draws.sens[d,])
    
  democrat.list[d] <- mean(pred.democrat.list)
  democrat.sens[d] <- mean(pred.democrat.sens)
  democrat.diff[d] <- democrat.list[d] - democrat.sens[d]
  
}

draws.list <- mvrnorm(n = n.draws, mu = par.list[1:k], Sigma = vcov.list[1:k,1:k])
draws.sens <- mvrnorm(n = n.draws, mu = par.sens, Sigma = vcov.sens)

for (d in 1:n.draws) {
  
  par.g <- draws.list[d, ]
    
  pred.republican.list <- logistic(X.republican.list %*% par.g)

  pred.republican.sens <- logistic(X.republican.sens %*% draws.sens[d,])
    
  republican.list[d] <- mean(pred.republican.list)
  republican.sens[d] <- mean(pred.republican.sens)
  republican.diff[d] <- republican.list[d] - republican.sens[d]

  
}

draws.list <- mvrnorm(n = n.draws, mu = par.list[1:k], Sigma = vcov.list[1:k,1:k])
draws.sens <- mvrnorm(n = n.draws, mu = par.sens, Sigma = vcov.sens)

for (d in 1:n.draws) {
  
  par.g <- draws.list[d, ]
    
  pred.independent.list <- logistic(X.independent.list %*% par.g)

  pred.independent.sens <- logistic(X.independent.sens %*% draws.sens[d,])
    
  independent.list[d] <- mean(pred.independent.list)
  independent.sens[d] <- mean(pred.independent.sens)
  independent.diff[d] <- independent.list[d] - independent.sens[d]

  
}

democrat.list.mean <- mean(democrat.list)
democrat.list.se <- sd(democrat.list)
democrat.sens.mean <- mean(democrat.sens)
democrat.sens.se <- sd(democrat.sens)
democrat.diff.mean <- mean(democrat.diff)
democrat.diff.se <- sd(democrat.diff)

republican.list.mean <- mean(republican.list)
republican.list.se <- sd(republican.list)
republican.sens.mean <- mean(republican.sens)
republican.sens.se <- sd(republican.sens)
republican.diff.mean <- mean(republican.diff)
republican.diff.se <- sd(republican.diff)

independent.list.mean <- mean(independent.list)
independent.list.se <- sd(independent.list)
independent.sens.mean <- mean(independent.sens)
independent.sens.se <- sd(independent.sens)
independent.diff.mean <- mean(independent.diff)
independent.diff.se <- sd(independent.diff)

democrat.list <- republican.list <- independent.list <-
  democrat.sens <- republican.sens <- independent.sens <-
  democrat.diff <- republican.diff <- independent.diff <- list()

for (college.list in 1:2) {

  college.set <- college.list - 1
  
  X.democrat.list <- as.matrix(subset(X.list, democrat == 1 & college == college.set))
  X.republican.list <- as.matrix(subset(X.list, republican == 1 & college == college.set))
  X.independent.list <- as.matrix(subset(X.list, democrat == 0 & republican == 0
                                         & college == college.set))
  
  X.democrat.sens <- as.matrix(subset(X.sens, democrat == 1 & college == college.set))
  X.republican.sens <- as.matrix(subset(X.sens, republican == 1 & college == college.set))
  X.independent.sens <- as.matrix(subset(X.sens, democrat == 0 & republican == 0
                                         & college == college.set))
  
  democrat.list[[college.list]] <- republican.list[[college.list]] <-
    independent.list[[college.list]] <-
    democrat.sens[[college.list]] <- republican.sens[[college.list]] <-
      independent.sens[[college.list]] <-
      democrat.diff[[college.list]] <- republican.diff[[college.list]] <-
        independent.diff[[college.list]] <-
        rep(NA, n.draws)
  
  draws.list <- mvrnorm(n = n.draws, mu = par.list[1:k], Sigma = vcov.list[1:k,1:k])
  draws.sens <- mvrnorm(n = n.draws, mu = par.sens, Sigma = vcov.sens)
  
  for (d in 1:n.draws) {
    
    par.g <- draws.list[d, ]
    
    pred.democrat.list <- logistic(X.democrat.list %*% par.g)
    
    pred.democrat.sens <- logistic(X.democrat.sens %*% draws.sens[d,])
    
    democrat.list[[college.list]][d] <- mean(pred.democrat.list)
    democrat.sens[[college.list]][d] <- mean(pred.democrat.sens)
    democrat.diff[[college.list]][d] <-
      democrat.list[[college.list]][d] - democrat.sens[[college.list]][d]
    
  }
  
  draws.list <- mvrnorm(n = n.draws, mu = par.list[1:k], Sigma = vcov.list[1:k,1:k])
  draws.sens <- mvrnorm(n = n.draws, mu = par.sens, Sigma = vcov.sens)
  
  for (d in 1:n.draws) {
    
    par.g <- draws.list[d, ]
    
    pred.republican.list <- logistic(X.republican.list %*% par.g)
    
    pred.republican.sens <- logistic(X.republican.sens %*% draws.sens[d,])
    
    republican.list[[college.list]][d] <- mean(pred.republican.list)
    republican.sens[[college.list]][d] <- mean(pred.republican.sens)
    republican.diff[[college.list]][d] <-
      republican.list[[college.list]][d] - republican.sens[[college.list]][d]
    
    
  }
  
  draws.list <- mvrnorm(n = n.draws, mu = par.list[1:k], Sigma = vcov.list[1:k,1:k])
  draws.sens <- mvrnorm(n = n.draws, mu = par.sens, Sigma = vcov.sens)
  
  for (d in 1:n.draws) {
    
    par.g <- draws.list[d, ]
    
    pred.independent.list <- logistic(X.independent.list %*% par.g)
    
    pred.independent.sens <- logistic(X.independent.sens %*% draws.sens[d,])
    
    independent.list[[college.list]][d] <- mean(pred.independent.list)
    independent.sens[[college.list]][d] <- mean(pred.independent.sens)
    independent.diff[[college.list]][d] <-
      independent.list[[college.list]][d] - independent.sens[[college.list]][d]
    
  }

}

democrat.diff.sim <-  democrat.diff[[1]] - democrat.diff[[2]]
republican.diff.sim <- republican.diff[[1]] - republican.diff[[2]] 
independent.diff.sim <-  independent.diff[[1]] - independent.diff[[2]]

educ.democrat.diff.mean <- mean(democrat.diff.sim)
educ.democrat.diff.se <- sd(democrat.diff.sim)
educ.republican.diff.mean <- mean(republican.diff.sim)
educ.republican.diff.se <- sd(republican.diff.sim)
educ.independent.diff.mean <- mean(independent.diff.sim)
educ.independent.diff.se <- sd(independent.diff.sim)
  
noncollege.democrat.list.mean <- mean(democrat.list[[1]])
noncollege.democrat.list.se <- sd(democrat.list[[1]])
noncollege.democrat.sens.mean <- mean(democrat.sens[[1]])
noncollege.democrat.sens.se <- sd(democrat.sens[[1]])
noncollege.democrat.diff.mean <- mean(democrat.diff[[1]])
noncollege.democrat.diff.se <- sd(democrat.diff[[1]])

noncollege.republican.list.mean <- mean(republican.list[[1]])
noncollege.republican.list.se <- sd(republican.list[[1]])
noncollege.republican.sens.mean <- mean(republican.sens[[1]])
noncollege.republican.sens.se <- sd(republican.sens[[1]])
noncollege.republican.diff.mean <- mean(republican.diff[[1]])
noncollege.republican.diff.se <- sd(republican.diff[[1]])

noncollege.independent.list.mean <- mean(independent.list[[1]])
noncollege.independent.list.se <- sd(independent.list[[1]])
noncollege.independent.sens.mean <- mean(independent.sens[[1]])
noncollege.independent.sens.se <- sd(independent.sens[[1]])
noncollege.independent.diff.mean <- mean(independent.diff[[1]])
noncollege.independent.diff.se <- sd(independent.diff[[1]])

college.democrat.list.mean <- mean(democrat.list[[2]])
college.democrat.list.se <- sd(democrat.list[[2]])
college.democrat.sens.mean <- mean(democrat.sens[[2]])
college.democrat.sens.se <- sd(democrat.sens[[2]])
college.democrat.diff.mean <- mean(democrat.diff[[2]])
college.democrat.diff.se <- sd(democrat.diff[[2]])

college.republican.list.mean <- mean(republican.list[[2]])
college.republican.list.se <- sd(republican.list[[2]])
college.republican.sens.mean <- mean(republican.sens[[2]])
college.republican.sens.se <- sd(republican.sens[[2]])
college.republican.diff.mean <- mean(republican.diff[[2]])
college.republican.diff.se <- sd(republican.diff[[2]])

college.independent.list.mean <- mean(independent.list[[2]])
college.independent.list.se <- sd(independent.list[[2]])
college.independent.sens.mean <- mean(independent.sens[[2]])
college.independent.sens.se <- sd(independent.sens[[2]])
college.independent.diff.mean <- mean(independent.diff[[2]])
college.independent.diff.se <- sd(independent.diff[[2]])

## generate figure

pdf("figure3.pdf", width = 6.5, height = 4)
par(mfrow = c(1, 2), cex = 0.7, mar = c(3, 4, 1, 1), las = 1)
plot(c(0.1, 0.2, 0.3), c(democrat.list.mean, democrat.sens.mean, democrat.diff.mean), pch = 16, xlim = c(0, 1.2), ylim = c(-0.5, 1.2), 
     xlab = "", ylab = "Estimated Proportion / Difference in Proportions", main = "", axes = FALSE)
points(c(0.5, 0.6, 0.7), c(republican.list.mean, republican.sens.mean, republican.diff.mean), pch = 17)
points(c(0.9, 1.0, 1.1), c(independent.list.mean, independent.sens.mean, independent.diff.mean), pch = 15)
axis(side = 2, at = seq(from = -0.4, to = 1, by = 0.2))
axis(side = 1, at = c(0.2, 0.6, 1.0), tick = FALSE,
     labels = c("Democrats", "Republicans", "Independents"), cex.axis = 0.8)

critical <- abs(qnorm(0.025))
lines(c(0.1, 0.1), c(democrat.list.mean - critical*democrat.list.se, democrat.list.mean + critical*democrat.list.se)) 
lines(c(0.2, 0.2), c(democrat.sens.mean - critical*democrat.sens.se, democrat.sens.mean + critical*democrat.sens.se)) 
lines(c(0.3, 0.3), c(democrat.diff.mean - critical*democrat.diff.se, democrat.diff.mean + critical*democrat.diff.se))

lines(c(0.5, 0.5), c(republican.list.mean - critical*republican.list.se, republican.list.mean + critical*republican.list.se)) 
lines(c(0.6, 0.6), c(republican.sens.mean - critical*republican.sens.se, republican.sens.mean + critical*republican.sens.se)) 
lines(c(0.7, 0.7), c(republican.diff.mean - critical*republican.diff.se, republican.diff.mean + critical*republican.diff.se))

lines(c(0.9, 0.9), c(independent.list.mean - critical*independent.list.se, independent.list.mean + critical*independent.list.se)) 
lines(c(1.0, 1.0), c(independent.sens.mean - critical*independent.sens.se, independent.sens.mean + critical*independent.sens.se)) 
lines(c(1.1, 1.1), c(independent.diff.mean - critical*independent.diff.se, independent.diff.mean + critical*independent.diff.se))

abline(h = 0, lty = "dashed")

text(0.1, 0.88, "List", cex = 0.8)
text(0.2, 0.4, "Direct", cex = 0.8)
text(0.3, 0.73, "Difference\nList - Direct", cex = 0.8)

plot(c(0.1, 0.2, 0.3), c(college.democrat.diff.mean, noncollege.democrat.diff.mean, educ.democrat.diff.mean), pch = 16, xlim = c(0, 1.2), ylim = c(-0.5, 1.2), 
     xlab = "", ylab = "Estimated Difference in Proportions (social desirability bias)", main = "", axes = FALSE)
points(c(0.5, 0.6, 0.7), c(college.republican.diff.mean, noncollege.republican.diff.mean, educ.republican.diff.mean), pch = 17)
points(c(0.9, 1.0, 1.1), c(college.independent.diff.mean, noncollege.independent.diff.mean, educ.independent.diff.mean), pch = 15)
axis(side = 2, at = seq(from = -0.4, to = 1, by = 0.2))
axis(side = 1, at = c(0.2, 0.6, 1.0), tick = FALSE,
     labels = c("Democrats", "Republicans", "Independents"), cex.axis = 0.8)

critical <- abs(qnorm(0.025))
lines(c(0.1, 0.1), c(college.democrat.diff.mean - critical*college.democrat.diff.se, college.democrat.diff.mean + critical*college.democrat.diff.se)) 
lines(c(0.2, 0.2), c(noncollege.democrat.diff.mean - critical*noncollege.democrat.diff.se, noncollege.democrat.diff.mean + critical*noncollege.democrat.diff.se)) 
lines(c(0.3, 0.3), c(educ.democrat.diff.mean - critical*educ.democrat.diff.se, educ.democrat.diff.mean + critical*educ.democrat.diff.se))

lines(c(0.5, 0.5), c(college.republican.diff.mean - critical*college.republican.diff.se, college.republican.diff.mean + critical*college.republican.diff.se)) 
lines(c(0.6, 0.6), c(noncollege.republican.diff.mean - critical*noncollege.republican.diff.se, noncollege.republican.diff.mean + critical*noncollege.republican.diff.se)) 
lines(c(0.7, 0.7), c(educ.republican.diff.mean - critical*educ.republican.diff.se, educ.republican.diff.mean + critical*educ.democrat.diff.se))

lines(c(0.9, 0.9), c(college.independent.diff.mean - critical*college.independent.diff.se, college.independent.diff.mean + critical*college.independent.diff.se)) 
lines(c(1.0, 1.0), c(noncollege.independent.diff.mean - critical*noncollege.independent.diff.se, noncollege.independent.diff.mean + critical*noncollege.independent.diff.se)) 
lines(c(1.1, 1.1), c(educ.independent.diff.mean - critical*educ.independent.diff.se, educ.independent.diff.mean + critical*educ.democrat.diff.se))

abline(h = 0, lty = "dashed")

text(0.1, 0.55, "College", cex = 0.8)
text(0.2, 1.02, "Non-College", cex = 0.8)
text(0.3, 0.07, "Difference\nNon-Coll - Coll", cex = 0.8)

dev.off()


## ####
## Figure 4 - "Monte Carlo Evaluation of the Three Estimators"
## pg 22

## See the PA_Simulations.R file for simulation replication code.


## ####
## Table 6 - "Estimated Respondent Types for the List Experiment in the 1991 National Race and Politics Survey"
## pg 31
##
## Also replicates in-text analysis: conducts test with null hypothesis that there is no design effect
## pg 30

rm(list=ls())
require(list)

data(affirm)
data(race)

ict.test(affirm$y, affirm$treat, J = 3, gms = TRUE, pi.table = TRUE)

ict.test(race$y, race$treat, J = 3, gms = TRUE, pi.table = TRUE)

## ####
## Table 7 - "Estimated Coefficients from the Logistic Regression Models"
## pg 31

rm(list=ls())
require(list)

# Fit standard design ML model
# Replicates Columns 1-2

noboundary.results <- ictreg(y ~ age + college + male + south, treat = "treat",
		      	     J = 3, data = affirm, method = "ml", 
			     overdispersed = FALSE)

summary(noboundary.results)

# Fit standard design ML model with ceiling effects alone
# Replicates Columns 3-4

ceiling.results <- ictreg(y ~ age + college + male + south, treat = "treat", 
		   	  J = 3, data = affirm, method = "ml", fit.start = "nls",
			  ceiling = TRUE, ceiling.fit = "bayesglm",
			  ceiling.formula = ~ age + college + male + south)

summary(ceiling.results, boundary.proportions = T)

# Fit standard design ML model with floor effects alone
# Replicates Columns 5-6

floor.results <- ictreg(y ~ age + college + male + south, treat = "treat", 
	      	 	J = 3, data = affirm, method = "ml", fit.start = "glm", 
			floor = TRUE, floor.fit = "bayesglm",
			floor.formula = ~ age + college + male + south)

summary(floor.results, boundary.proportions = T)

# Fit standard design ML model with floor and ceiling effects
# Replicates Columns 7-8 

both.results <- ictreg(y ~ age + college + male + south, treat = "treat", 
	     	       J = 3, data = affirm, method = "ml", 
		       floor = TRUE, ceiling = TRUE, 
		       floor.fit = "bayesglm", ceiling.fit = "bayesglm",
		       floor.formula = ~ age + college + male + south,
		       ceiling.formula = ~ age + college + male + south)

summary(both.results, boundary.proportions = T)

## ####
## Figure 5 - "Statistical Power of the Proposed Test to Detect Design Effects"
## pg 34

## See the PA_Simulations.R file for simulation replication code.

##eof
